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ABSTRACT 

A^-body simulations predict that CDM halo-assembly occurs in two phases: 1) a fast 
accretion phase with a rapidly deepening potential well; and 2) a slow accretion phase 
characterised by a gentle addition of mass to the outer halo with little change in the 
inner potential well. We demonstrate, using one-dimensional simulations, that this 
two-phase accretion leads to CDM halos of the NFW form and provides physical in- 
sight into the properties of the mass accretion history that influence the final profile. 
Assuming that the velocities of CDM particles are effectively isotropised by fluctu- 
ations in the gravitational potential during the fast accretion phase, we show that 
gravitational collapse in this phase leads to an inner profile p(r) oc r~^. Slow accre- 
tion onto an established potential well leads to an outer profile with p{r) cx 7-^'^. The 
concentration of a halo is determined by the fraction of mass that is accreted during 
the fast accretion phase. Using an ensemble of realistic mass accretion histories, we 
show that the model predictions of the dependence of halo concentration on halo for- 
mation time, and hence the dependence of halo concentration on halo mass, and the 
distribution of halo concentrations all match those found in cosmological A^-body sim- 
ulations. Using a simple analytic model that captures much of the important physics 
we show that the inner r^^ profile of CDM halos is a natural result of hierarchical 
mass assembly with a initial phase of rapid accretion. 
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1 INTRODUCTION 

In the cold dark matter (CDM) paradigm of structure forma- 
tion, most of the cosmic mass is locked in virialised clumps 
called dark matter halos. Luminous objects, galaxies and 
clusters of galaxies, are assumed to form in the potential 
wells of these dark matter halos. High resolution A''-body 
simulations have shown that the density profiles of CDM 
halos can be fairly well described by a universal form. 
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where is a characteristic radius and ps is a characteris- 
tic density (Navarro, Frenk & White 1996, 1997; hereafter 
(NFW). The value of is often given in units of the virial 
radius and one over that value is referred to as the halo con- 
centration. There is still uncertainty about the exact value 
of the inner slope. While some simulations indicate that the 
inner logarithmic slope may be steeper than the NFW value, 
-1 (e.g. Moore et al. 1999, Ghigna et al. 2000; Fukushige & 
Makino 1997, 2001, 2003), others give slopes shallower than 
— 1 (Subramanian, Cen & ©striker 2000; Taylor & Navarro 
2001; Ricotti 2003). Jing & Suto (2000) found that CDM 
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halo profiles in their simulations do not have a completely 
universal form, with the inner slope changing from system 
to system. 

Until now, there has not been a natural explanation 
for the approximate 'universal' profile resulting from the 
gravitational collapse of a CDM density field. The pioneer- 
ing work by Gunn & Gott (1972) considered the collapse 
of uniform spherical perturbations of coUisionless cold dark 
matter in an expanding background. This simple model ex- 
plains some properties of virialised objects, such as their 
mean density and size, but does not describe the density 
profile of a collapsed object. Subsequently, the spherical col- 
lapse model was extended to incorporate realistic initial per- 
turbations (e.g. Gott 1975; Gunn 1977; Fillmore & Goldre- 
ich 1984; Bertschinger 1985; Hoffman & Shaham 1985; Ry- 
den & Gunn 1987; Ryden 1988; Zaroubi & Hoffman 1993). 
Bertschinger (1985) and Fillmore & Goldreich (1984) found 
that halos have similar asymptotic inner density profiles 
with p(r) oc and 7 ~ 2, for a wide range of initial 
density perturbations. Unfortunately, the inner slope pre- 
dicted by such models is much steeper than that found in 3- 
dimensional cosmological simulations. A number of authors 
(e.g. White & Zaritsky 1992; Ryden 1993; Sikivie et al. 1997; 
Subramanian 2000; Subramanian et al. 2000; Hiotelis 2002; 
Le Delliou & Henriksen 2003; Shapiro et al. 2004; Barnes et 
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al. 2005) noticed that tangential motions of particles may 
cause flattening of the inner profiles. A'^-body simulations 
by Huss et al. (1999) and Hansen & Moore (2004) showed 
that the inner density profile of a halo is correlated with the 
degree of velocity anisotropy. However, none of these has 
provided clear dynamical mechanism for the origin of the 
'univeral' p(r) oc inner profile observed in cosmological 
A''-body simulations. 

Lynden-Bell (1967) proposed that, as long as the ini- 
tial condition for the collapse is clumpy, a final equilibrium 
state with a universal profile may be achieved as a result 
of violent relaxation that may cause a complete mixing of 
particles in phase space. Numerical simulations have shown 
that such collapses indeed produce a universal inner profile 
(van Albada 1982), but the relaxation is incomplete, in the 
sense that there is still a significant correlation between the 
final and initial states of a particle. Tremaine et al. (1986) 
provide statistical constraints on equilibria resulting from 
violent relaxation. In the cosmic density field predicted by 
a CDM model, the perturbations responsible for the for- 
mation of dark matter halos are expected to be clumpy, 
and so violent relaxation is expected to play some role in 
the formation of dark matter halos. However, as shown by 
A'^-body simulations, the density profiles of virialised dark 
matter halos do depend on the formation histories of dark 
halos (e.g. NFW; Klypin et al. 2001; Bullock et al. 2001; 
Eke et al. 2001; Wechsler et al. 2002; Zhao et al. 2003a,b; 
Tasitsiomi et al. 2004; Diemand et al. 2005). NFW conjec- 
tured that the dependence of halo concentration parameter 
on halo mass could be explained by the assembly of more 
massive halos at later times than lower-mass halos. Later, 
Wechsler et al. (2002) and Zhao et al. (2003a;b) found that 
the concentration of a halo depends on the assembly history. 
Together, these results suggest that CDM initial conditions 
play an important role in determining halo density profiles. 

In this paper, we explore the relation between the den- 
sity profile and the assembly histories of dark halos in CDM 
models. Our goal is to understand the physical processes 
that shape the density profiles of CDM halos, further moti- 
vated by the recent finding that the mass accretion histories 
of CDM halos show remarkable regularity. Based on high- 
resolution A'^-body simulations, the mass accretion history of 
a halo in general consists of two distinct phases: an early fast 
phase and a later slow phase (Wechsler et al. 2002; Zhao et 
al. 2003a,b; Li et al. 2005). As shown in Zhao et al. (2003a), 
the fast accretion phase is dominated by major mergers char- 
acterised by a rapid deepening of the halo potential. The 
slow accretion phase is characterised by only weak changes 
in the gravitational potential well. We will demonstrate that 
the 'universal' NFW profile results from gross properties of 
the mass accretion histories of dark matter halos. Using a 
simple spherical collapse model, we show that the inner pro- 
file is established in the fast accretion phase. Its key features 
are rapid collapse (in a small fraction of a Hubble time) with 
rapid changes in the potential that may isotropise the ve- 
locity field. The outer profile is dominated by particles that 
are accreted in the slow-accretion phase onto an existing 
central object. Our model incorporates these two phases of 
CDM accretion history to explain the appearance of 'univer- 
sal' density profile in A^-body simulations. In the absence of 
slow accretion, the outer profile would approach p oc r~*. 



Our findings suggest that mass accretion history plays a 
crucial role in structuring halos. 

The outline of this paper is as follows. In ^ we briefiy 
describe CDM halo mass accretion histories, our procedure 
for realizing initial conditions for the formation of dark mat- 
ter halos, and present our one-dimensional algorithm for 
simulating the collapse of dark matter halos. Our simula- 
tion results are described in JH] In SjJl ^6 use models with 
simple power-law initial perturbations to understand how 
NFW-like profiles are produced in CDM models. Finally, in 
we further discuss and summarise our results. 



2 INITIAL CONDITIONS AND METHODS 

For a given halo, the mass accretion history specifies how 
much mass is added to the halo as a function of time. Nu- 
merical simulations and analytical models of halo formation 
in CDM models reveal that the mass accretion histories of 
CDM halos are remarkably regular. Wechsler et al. (2002) 
found that the accretion history of a CDM halos from TV- 
body simulations can be described roughly by the following 
parametric form: 



M (a) = Mo exp 



(2) 



where a is the expansion scale factor of the universe, and Mo 
is the virial mass of the halo at a final time where a = ao. 
The formation history is characterised by a single parame- 
ter He. This characteristic scale factor a — is the point 
when the logarithmic mass accretion rate, d log M/d log a, 
falls below a critical value S = 2 (Wechsler et al. 2002). 
Similar results were found by Zhao et al. (2003a; b) using 
high-resolution simulations and by van den Bosch (2002) 
using extended Press-Schechter theory. 

In a cosmological spherical collapse model, the shell col- 
lapse time is determined by the mean initial over-density 
within the mass shell; the mass within a mass shell collapses 
when the average linear over-density (calculated using lin- 
ear perturbation theory) reaches Sc ~ 1.686. Therefore, for 
a given mass accretion history, we can then construct the 
corresponding initial density perturbation profile that re- 
produces this history. For a spherical perturbation of mass 
M that collapses at a redshift z, its linear over-density at 
the initial time Zi is given by 



5^{M) = 1.686 



D{z) 



(3) 



where D{z) is the linear growth factor. In our calculation, 
we use the fitting formula by Carroll et al. (1992), 
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(5) 



where ^m{z) and ^h{z) are the density parameters of non- 
relativistic matter and of the cosmological constant at red- 
shift z, respectively. At redshift Zi, the radius rt of the sphere 
that encloses mass M is 
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3M 

ATvp{zi)[l + S,{M)] 



1/3 



(6) 



where pi{zi) = Pcrit.of^Af (0)(1 + Zi)'^, with pcrit.o the criti- 
cal density of the universe at « = 0. Assuming that Zi is 
sufficiently large that the initial collapse is linear, equation 
^ relates the enclosed mass M{z) to its overdensity. Equa- 
tion JSJ then determines its radius at Zi. Choosing a mass 
partition A/i < M2 < ■ ■ ■ < Mm, one can then determine 
the radii ri < r2 < • • ■ < tn. We choose equal mass shells 
m = Mj+i — Mj for our simulations. Note that this proce- 
dure described by equations ^ and Q is independent of 
the mass scale Mq in equation ^ and, therefore, a single 
simulation may be scaled to any value of Mq ex post facto. 

We use one-dimensional particle simulations to explore 
the dynamical evolution of the collapse for a given mass ac- 
cretion history. Note that even though the simulations are 
one-dimensional we can include the effects of angular mo- 
mentum (see eq.|5] below). With a correction for the cosmo- 
logical constant, we approximate the gravitational accelera- 
tion of the fcth mass shell (particle) by 



(r2+Q2)3/2 



Mt = 



E 



(7) 
(8) 



where Ho is Hubble's constant at the present time, JIa is 
the density parameter of the cosmological constant, is 
the radius of the shell, and a is a softening length. In every 
simulation, we use a softening length a — 0.0005r„, where 
r„ is the virial radius of the halo at the present time. This 
scale is much smaller than any scale of interest. We assign 
each shell a specific angular momentum Jk; Jk = defines 
pure radial infall. The effective acceleration including the 
centrifugal force is then 



ak = gk -\ 3 



(9) 



Since in spherical symmetry the gravitational force on 
a mass shell is determined by the mass enclosed by the mass 
shell, we calculate the force by sorting particles according to 
their radii. Owing to the finite number of shells, when two 
shells cross they experience a gravitational force discontinu- 
ity since then the enclosed mass instantly changes. To reduce 
such effects, we introduce another type of force softening. 
Instead of assuming that each shell is infinitesimally thin, 
we assume that it has a constant density with finite thick- 
ness, which we choose to be the distance between the inte- 
rior and exterior neighbouring shells. Using such a softening 
method, crossing shells gradually change their enclosed mass 
and hence the gravitational forces change smoothly. 

We use a time-symmetric symplectic leapfrog integrator 
(Quinn et al. 1997; Springel 2005) to solve the equations of 
motion. At each time step, we calculate the dynamical time 
for every shell and we choose the next time step to be smaller 
than the shortest dynamical time of all the shells, i.e. 



AW=mm<jc.^ 2GM. 



5rl 



(10) 




where Cd is control parameter, which we set to be 0.001. Each 
simulation starts from an initial condition, which specifies 



Figure 1. The diamonds are the binned density profile of a simu- 
lated dark matter halo at 2 = with Poisson error bars. Particles 
are assumed to have pure radial motion in this simulation. Note 
that the outer density profile has p oc r^'^, while the inner profile 
has p (X r~^. The dashed curve shows a NFW profile for compar- 
ison. 



the position and velocity of each particle at a chosen high 
redshift. The position for each particle at the initial red- 
shift Zi follows from the mass accretion history as described 
above. The initial velocity consists of two components: the 
Hubble expansion Vi{M) — ri{M)H{zi) and the peculiar ve- 
locity. We use linear perturbation theory to relate the initial 
peculiar velocity of a mass shell to Si{M). We start our sim- 
ulations at Zi = 200, early enough for linear theory to be 
valid. 

All our simulations assume Hm ~ 0.3 and Qa = 0.7. We 
use 10^ equal mass particles to simulate the formation of a 
single dark matter halo, and we have tested that this number 
is sufficiently large to achieve numerical convergence over 
the scales in which we are interested. In each simulation the 
total mass of all the particles is 1.2 times the final virial mass 
of the halo at z — 0. The extra mass maintains an ambient 
environment to follow the late stages of halo formation. As 
a test of our code we have reproduced the self-similar results 
of both Fillmore & Goldreich (1984) and Nusser (2001). We 
varied the particle number by a factor of 10, both larger and 
smaller, and could still satisfactorily reproduce these results. 



3 RESULTS 

We first consider a model with ac — 0.4 corresponding to 
a formation redshift of ~ 1.5 with pure radial motion 
{Jk = in eq. |5J. The density profile of the halo at z = 
is shown in Figure Q The final z = density profile has an 
inner logarithmic slope of —2 and an outer slope of —3. The 
predicted inner profile is much steeper than a NFW pro- 
file, although the model matches a NFW profile in the outer 
parts. We find that purely radial collapse simulations with 
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varying values of ac all produce halo profiles with an inner 
logarithmic slopes of —2 and outer slopes of —3. Differing 
values of ac affect only the transition radius between these 
two slopes; larger values of Qc (later formation times) leads 
to a larger transition radius, i.e. a lower concentration. Ra- 
dial infall model alone cannot reproduce a NFW profile in 
the inner parts because it does not accurately represent the 
dynamics of fast accretion. The early fast accretion phase of 
a halo is dominated by major mergers and the depth of the 
potential well associated with the main progenitor deepens 
rapidly with time (Zhao et al. 2003a). Frequent scattering 
by potential fiuctuations in this phase is expected to effec- 
tively isotropise the orbits. As a result, dark matter particles 
will acquire a significant amount of angular motion as they 
are accreted by the halo, which is not included in the purely 
radial calculation. However, the radial infall model does ap- 
pear to successfully describe accretion from large distances 
onto an existing central mass in the slow accretion regime of 
halo formation and we will see in 21 that this does explain 
the agreement of the model with a NFW density profile in 
the outer parts. To simulate the fast accretion process more 
accurately using our one-dimensional approach, we consider 
a model that includes the effect of isotropic velocity disper- 
sions. We trace the motion of each particle assuming pure 
radial motion into a radius of i?t, and at this point, we ran- 
domly assign a tangential component of velocity to the parti- 
cle, keeping the kinetic energy unchanged. As the simulation 
evolves further, the angular momentum of each particle is 
conserved after an 'isotropisation' event. We choose Rt to 
be one half of the turn-around radius of each particle. Al- 
though this prescription seems somewhat arbitrary, a differ- 
ent choice only moderately shifts the radial scale. In detail, 
we incorporate this into our isotropisation prescription as 
follows. For a given mass shell M, we assume the ratio be- 
tween the radial velocity dispersion, Gt, and the tangential 
velocity dispersion at, has the form. 



(11) 



where ra is the characteristic scale of the halo demarcat- 
ing fast and slow accretion and 13 controls the shape of the 
anisotropy profile. The radial and tangential velocities are 
randomly sampled from Gaussian distributions. The ratio 
of these two random numbers are then used to partition the 
kinetic energy of the particle into radial and tangential com- 
ponents. For mass shells with Rt <^ ra, the velocity distribu- 
tion is nearly isotropic, while for those with Rt 3> ra radial 
orbits dominate. Note that a larger /3 produces a sharper 
transition between radial and isotropic orbits. We take P = 2 
as our fiducial model and will describe the effect of changing 
the value of /3 below. 

Since Rt is roughly the virial radius of a mass shell, we 
take ra = rv{ac), where r„(ac) is the virial radius of the halo 
at the transition time between the fast accretion phase and 
the slow accretion phase. Figure |21 shows the final halo den- 
sity profile of such a simulation. The resulting inner density 
profile is much flatter than the p oc r~^ we find for pure 
radial motions, and now over a large range of radius it can 
be well fit by a NFW profile. The inner profile is dominated 
by particles accreted in the fast accretion regime, and the 
shallower profile owes to the non-radial motions of these par- 
ticles. These results are not particular to the mass accretion 
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Figure 2. The density profile of a dark matter halo (diamonds) in 
a simulation where particles accreted in the fast accretion phase 
are assumed to have isotropic velocity dispersion (see text for 
details). The simulated density profile can be fit by a NFW profile 
(dashed curve). The dotted curve is the result of the best fit to 
the data by equation 1121 . and the best- fit value for 7 is 1.21. 
Squares show the density profile for particles that are accreted in 
the fast accretion phase, while triangles show the density profile 
for particles accreted in the slow accretion phase. Note that the 
outer part of the halo profile is dominated by particles in slow 
accretion, while the inner profile is dominated by particles in the 
fast accretion phase. 



history used in this example, and indeed our simulations 
with other mass accretion histories (see below) all lead to 
similar results. 

We check the sensitivity to our velocity structure as- 
sumptions by varying the values of the two parameters 
and f3 in equation IIH . We first consider two alternatives to 
the fiducial model — r„(ac): ra = 0.5r„(ac) and 2r„(ac) 
with fixed [3 = 2. The density profiles given by these two 
models are shown in the left panel of Figure 01 The den- 
sity profile changes very little even though the value of ra 
changes by a factor of 4. Next, we hold ra fixed at our fidu- 
cial value of r-„(ac) and vary the value of /3 from 0.5 to 10. 
The resulting density profiles are shown in the right panel 
of Figure 13 The final profile is also insensitive to changes in 
/3. 

In summary, if particles accreted in the fast accretion 
regime are assumed to have an isotropic velocity dispersion, 
we find that the inner slope of the halo density profile always 
lies in the range between —1 and —1.2. Apparently, the in- 
ner p (X r^^ profile is a generic result of the fast collapse 
phase and the isotropic velocity field generated during such 
a collapse. 

The results obtained above can be compared to those 
obtained earlier by Avila- Reese et al. (1998) and Ascasibar 
et al. (2004). These authors examined the density profiles of 
dark matter halos produced by spherical collapse in a CDM 
density field, assuming that the orbits retain a constant el- 
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Figure 3. The density profiles for models with different choices of and f) [see equation 1111 for definition]. In the left panel, results 
are shown for models where /3 is fixed to be 2 but ra changes from rv{ac) (diamonds) to 0.5rv{ac) (squares) and 2rv{ac) (triangles). In 
the right panel, ra is fixed to be ry(ac) but the value of (3 changes from 2 (diamonds) to 0.5 (squares) and 10 (triangles). The solid curve 
shows a NFW profile, while the dashed curve shows the profile (ean. IT2I with 7 = 1.2. 



lipticity. These models can produce NFW-like profiles pro- 
vided that the constant ellipticity is properly chosen. How- 
ever, as we describe in ^ spherical collapse with constant 
ellipticity orbits can produce a wide range of inner profiles 
depending on initial conditions. A good match between their 
models and the NFW profile requires a fine tuning of initial 
conditions. Shapiro et al. (2004) also explored the impor- 
tance of CDM accretion history using one-dimensional sim- 
ulations and very similar arguments to ours. Unfortunately, 
they use a fluid approach that solves Jeans-like moments of 
the coUisionless Boltzmann equation and, presumably, the 
elimination of any possible asymmetric velocity distribution 
prevented them from finding our p oc result. 

In contrast to these models, our consideration is based 
on realistic CDM halo formation histories. We demonstrate 
that, for all such formation histories, the early fast accretion 
of dark matter that may effectively generate an isotropic 
velocity dispersion leads naturally to p oc in the in- 
ner parts. Our isotropisation assumption is supported by 
the measurement of halo velocity dispersions in cosmologi- 
cal iV-body simulations, which becomes progressively more 
isotropic towards the inner part of dark halos (e.g. Eke 
et al. 1998; Colin et al. 2000; Fukushige & Makino 2001; 
Hansen & Moore 2004). 

We can study the dependence of halo structure on 
the 'formation time', characterised by Uc, using our one- 
dimensional simulations. We fit the resulting density profiles 
both with a NFW profile and one with the following more 
general form: 
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Figure 4. Halo concentration as a function of ac, which char- 
acterises the halo formation time. Triangles show the results 
where simulated halo profiles are fit with a NFW profile and 
the diamonds are the results when fit with equation 1121 . The 
solid curve shows the prediction of a model proposed by Zhao et 
al. (2003a) based on cosmological A'^-body simulations, while the 
dashed curve shows the model prediction given by Wechsler et 
al. (2002). 



p{r) = 



(r/rs)'' [l + {r/rs)f-^'> ' 
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where 7 is the inner slope of the profile. Figure 21 shows 
the best fit value of the concentration parameter, defined 
as c = r^/rs, as a function of flc. As one can see, halos 
with lower Oc, i.e. earlier formation times, are more con- 
centrated. However, for halos that are still in the fast ac- 
cretion phase, i.e. with Qc > 1, the concentration is inde- 
pendent of Sc. The solid curve shows the concentration- 
formation time relation obtained by Zhao et al. (2003a;b) 
using high-resolution A'^-body simulations. Our results as- 
suming a NFW profile match this relation extremely well. 
The dashed line shows the model proposed by Wechsler et 
al. (2002). This model agrees well with our results for ac < 1, 
but underestimates the concentration for halos with Uc > 1. 
We refer readers to Zhao et al. (2003a;b) for a detailed 
discussion about the discrepancy between their results and 
those obtained by Wechsler et al. (2002). Since it is known 
that the concentration-formation time relation is the ori- 
gin of the mass-concentration relation (Wechsler et al. 2002, 
Zhao et al. 2003a), our results will also reproduce the mass- 
concentration relation (at a given redshift) as well as the 
concentration-redshift relation (at fixed mass) found in the 
cosmological A-body simulations. 

So far, our discussion is based on simulations that as- 
sume the smooth accretion history given by equation 
Although this smoothed form is a good description of the 
accretion history averaged over an ensemble of simulated 
dark halos, an individual accretion history exhibits details 
that are not described by equation Q. In some cases, equa- 
tion @ even fails to describe the overall shape of the mass 
accretion history. We show such an example in the left panel 
of Figure |5] A fit that emphasises the accretion history be- 
fore the big jump at a = 0.45, shown by the short dashed 
curve, differs from a fit that emphasises the history after 
that time. To see how such differences affect the density 
profiles, we carried out one-dimensional simulations using 
realistic halo accretion histories generated by PINOCCHIO, 
a Lagrangian code developed by Monaco et al. (2002). The 
statistical properties of the mass accretion histories gener- 
ated using this code agrees with those obtained using cos- 
mological A-body simulations (Li et al. 2005) . The resulting 
simulated halo density profiles can all be well described by 
the NFW form. For many cases where equation (I^J is a good 
fit to the overall accretion history, the simulated density pro- 
files using real accretion histories are very similar to those 
using the corresponding fits with equation ||5J . Even for cases 
where the accretion history is poorly described by equation 

as the one shown in Figure the final halo profile 
can still be well fit using a NFW form (see Fig.jSJ)). In such 
cases, however, the formation time defined using equation 
J^J is uncertain. 

The difference in the detailed mass accretion history in- 
troduces a scatter in the distribution of the concentration 
parameter c (Wechsler et al. 2002). If the density profile is 
determined by its mass accretion history, then we should 
be able to reproduce the distribution of c using a random 
sample of realistic mass accretion histories. To perform this 
experiment, we have randomly chosen 50 mass accretion his- 
tories generated with PINOCCHIO for halos with masses in 
the range mass accretion histories are 

used to generate initial conditions for our one-dimensional 
simulations. Figure |S] shows the distribution of the concen- 
tration parameter, obtained by fitting the simulated profiles 
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Figure 6. The histogram shows the distribution of halo concen- 
tration obtained from an ensemble of 50 randomly chosen mass 
accretion histories for halos with masses in the range 10^-*^ Mq to 
10^^ Mq at the present time. The solid curve shows a log-normal 
distribution with a median equal to 15.0 and a dispersion (in 
logc) equal to 0.12. 



with the NFW form. The distribution can be roughly de- 
scribed by a log-normal, in agreement with the cosmological 
A-body results (e.g. Jing 2000; Bullock et al. 2001). The dis- 
persion in log c, a = 0.12, is also very close to that found in 
cosmological A-body simulations (e.g. Wechsler et al. 2002), 
reinforcing our finding that the density profile of a halo is 
largely determined by its mass accretion history. 



4 WHAT DETERMINES THE DENSITY 
PROFILES OF DARK MATTER HALOS? 

In the last section, we have shown that many of the struc- 
tural properties of CDM halos found in cosmological A-body 
simulations can be understood in terms of halo mass accre- 
tion histories. However, why does gravitational collapse with 
such initial conditions always produce halo profiles that fol- 
low a universal form? Is the universal profile a result of the 
fact that the initial conditions represented by CDM mass 
accretion histories are special, or is it more generic in the 
sense that it can be produced by a wide range of initial 
conditions? 

To answer this question, we consider the collapse of 
generic initial perturbations with the form 



oc M(r) 



(13) 



SM{r) 
M{r) 

where e is a constant that controls the mass accretion rate. 
Because a mass shell collapses when 5M/M ~ 1.68, the 
mass accretion history implied by equation 11311 is M{z) oc 
D^^^{z), where D{z) is the linear growth factor. Since the 
circular velocity of a halo Vc is related to its mass, M ~ 
Vc /H{z), where H{z) is the Hubble constant at redshift z, 



© 2005 RAS, MNRAS OOO.ITHTTI 



Profiles of dark halos 7 



O'.C'l 





1 1 1 1 1 


















/ / 


y 
















/ /T 






1 /f 






1 \ 




r / 
/ 

/ 

/ 

/ 

/ 


Ji 

III 
jji 

1 , , , 1 


, , , 1 . , , 1 , . , 




0.0 0.2 0.4 0.6 0.8 

a 



1.0 



0.01 0.10 1.00 

r /r^ 



Figure 5. The diamonds in panel (a) show a mass accretion history that cannot be well fit using equation J^J- The long (short) dashed 
curve is the result of fit that emphasises the recent (past) history at a > 0.45 (a < 0.45). The solid curve represents a compromise 
between the recent and past histories. The halo profiles corresponding to the actual mass accretion history and these three possible fits 
are shown in panel (b). 



we can write Vc{z) oc H^^'^{z)D^^^'^{z). For simplicity, we 
consider an Einstein-de Sitter universe. In this case, D oc 
H~^^^ , and therefore 



M xH 



-2/3e 



VcxH 



(l-2/3e)/3 



(14) 



Note that Vi? is a measure of the depth of potential well as- 
sociated with the halo. For e = 1/6, Vc oc {M oc H"^). 
Therefore, e = 1/6 separates the isotropisation regime from 
calm accretion, i.e. Vc changes by an order of unity or more 
in a Hubble time for e < 1/6. For e = 2/3, M oc H'^ 
(Vc — constant) and, similarly, e — 2/3 separates fast ac- 
cretion from slow accretion in terms of the mass accretion 
rate. For e — > 0, perturbations on different mass scales have 
the same amplitude, and all mass scales collapse simultane- 
ously, leading to fast increases in both mass and potential 
well depth. In contrast, for e ^ 1 the perturbation ampli- 
tude declines rapidly with mass scale, hence M increases lit- 
tle while the potential well decays as the universe expands. 
This allows us to study cases with vastly different collapse 
histories by changing the value of e over a large range. 

In an Einstein-de Sitter universe where the expansion 
factor a is a power law of time, the collapse of perturba- 
tions described by equation 11311 admits similarity solutions. 
Assuming spherical symmetry and pure radial motion, Fill- 
more & Goldreich (1984) show that the collapse develops an 
asymptotic density profile p oc r ' in the inner region, with 
7 = 2 for < e < 2/3, and 7 = 9e/(H-3e) for e > 2/3. These 
solutions are shown in Figure |7| as the long-dashed curve. 
Self-similar models can also be constructed for non-radial 
motions (Nusser 2001). The specific angular momentum of 
a particle can be specified as 




0.5 



0.0 



isotropic 

K — 



constant / 



0.01 0.1 



10 100 



Figure 7. The relation between the inner logarithmic slope, 7 
(p oc r~^), as a function of the exponent of the initial pertur- 
bation defined in equation 1131 . The long dashed curve shows 
the solution of the model with pure radial motion and the dot- 
ted curve shows the solution when all particles have the same 
jr. The solution with an isotropic velocity dispersion is shown as 
the solid curve, and is compared with the result obtained directly 
from numerical calculations (crosses). All the three curves overlap 
for e > 2/3. 



J = J'VGA/tarta 



(15) 



© 2005 RAS, MNRAS OOO.miTTl 



8 



where rta is the turnaround radius of the particle, and 
Mta is the mass interior to rta- If J is the same con- 
stant for all particles, the problem admits similarity solu- 
tions and the asymptotic inner slope of the density profile 
is 7 = 9e/(l + 3e) for all e > (Nusser 2001). This solution 
is shown as the dotted curve in Figure □ For e > 2/3, this 
solution is the same as that for the model with pure radial 
infall. The asymptotic value of 7 is as e ^ 0. The depen- 
dence of 7 on e for a constant J is steep near 7 ~ 1. Thus, 
to produce a NFW inner slope in this model requires tuned 
initial conditions. 

For our models in the velocity dispersion of particles 
in the early collapse phase is isotropic. Thus, the quantity 
J has the following distribution at fixed energy: 

J 



P{J) 



(16) 



For the slow accretion phase, we expect the results of pure 
radial motion to obtain. Thus, for e > 2/3, the asymp- 
totic slope is 7 = 9e/(l + 3e). For e < 2/3, we can also 
use a simple model to understand the results obtained in 
the last section. Consider all particles with a given . In 
general, we can write the final density profile of these par- 
ticles as Apj{r) oc {AMj /rj)F{r/r j), where AMj is 
the total mass of particles with J in the range J ± AJ, 
and rj- is a characteristic scale in the density profile. The 
quantity AMj/rj is the density scale. Since the only scale 
in the problem is the current turnaround radius rout, and 
since non-radial motion is expected to be important only 
for r < J rout, we expect rj oc J rout- Thus, 



Apj{r) oc 



M P{J)AJ _ 

' out 



\Jr out J 



(17) 



If we neglect the interaction between mass shells of different 
J , the total density can be written as 



(18) 



M^.^i:^..>-.«i^/^.(^) 

On scales where the effect of angular momentum is im- 
portant, i.e. r ^ Jrout, the density profile is expected to 
be F(r) oc r~'^^^''-^^^'-^ . Conversely, for r 3> Jr out the parti- 
cles are in nearly radial motion and one expects F(r) oc 
(for e < 2/3). Hence, we may approximate the form of F as 
F{x) = I jlx^ -I- a;®'/f^+^^']. Inserting this and equation llTCl 
into equation H18|l . we obtain 



p(t) « ^'^—Qi^-.Toutlr) , 



where 



Q (e, rout/r) 



r/2 



(rout/r)de 



1 + [(rout/r) sine]- ' 



with 



2- 3e 



(19) 



(20) 



(21) 



l + 3e 

The inner density profile is given by the r dependence of 
Q for r/rout — > 0. The solid curve in Figure [7| shows the 
resulting 7-e relation. The characteristics of this relation can 
be understood as follows. For e < 1/6, then q > 1. The 
integration in Q is dominated by small 6, and so we can 
replace sinO by 9. The function Q is independent of r and 
p oc r"^. For 1/6 < e < 2/3, then < a < 1 and the 



integration is now dominated by sin 9 > r/rout - The function 
Q oc r""^ and, therefore, 7 — 9e/(l + 3e). This relation 
between 7 and e holds also for e > 2/3, as discussed above. 

To check the accuracy of the above simple model, we 
use a numerical calculation to solve for the inner density 
profile as a function of e. For any radius r, the total mass 
within it can be written in two parts. 



mrir) = mp{r) + mt{r) . 



(22) 



Here rrip is the mass of all particles with apocenter smaller 
than r. We call these particles 'permanent' contributors be- 
cause they always contribute to rriTir)- The mass rat in the 
above equation is the contribution of particles with apocen- 
ter ra larger than r but with a smaller pericenter r^,. These 
are 'temporary' contributors because they only spend part 
of their orbital times within r. Let P(r|rfe) be the fraction 
of time that a 'temporary' particle with turnaround radius 
rfc spends inside radius r. The total mass contributed by 
temporary particles can be written as 

mt{r) = / P(r|rfe)drn(rfe), (23) 

where the measure includes all orbits that pass through the 
surface at r. By definition. 



P(r\rk) = 



dr' 

,, Vk{r')i 



Ar' 

Vk(r') 



(24) 



where Vkir') is the radial velocity of a particle with turn- 
around radius at a radius r' . Note that Vk{r) as a function 
of r depends on the current mass profile mT(r), and so the 
mass profile has to be solved by iteration. The radial velocity 
can be written as 



Vk{r) = ^/2 [Ek - <E-(r) - Jl/{2r^)\ 



1/2 



(25) 



where and <& are, respectively, the total energy and grav- 
itational potential of the particle, and Jk is the specific an- 
gular momentum. The increase of mass within the apocenter 
of a particle can change the orbit of the particle. Thus, we 
must recompute such changes at each step of the iteration. 
Assuming that angular momentum is conserved, the apoc- 
enter ra of a particle after a given iteration is related to the 
apocenter r^ before by mT{ra)ra ~ m'j'{r-'a)r'a, where m(p 
and TUT are the mass profiles before and after the iteration 
step. The iteration starts from an initial condition where 
each particle is at its turnaround radius, and we denote the 
corresponding profile by Mta(rta). Since the turnaround ra- 
dius rta of a mass shell is related to its initial radius r^ by 
rta oc ri{SM/M)~^ in an Einstein-de Sitter universe and 
since M oc rf, one can show that 



Mta(rta) OC rra , u = 3/(3e + 1) . 



(26) 



We have applied the above numerical calculation to 
models with pure radial infall and models where J is the 
same for all mass shells. The results match those given by 
the self-similarity solutions. The crosses in Figure |7| show 
the 7-e results for this calculation using an isotropic veloc- 
ity dispersion given by equation II16II . The numerical results 
match the simple analytical model presented above remark- 
ably well. We have also applied our one-dimensional code 
to this model and the results are similar to those obtained 
here. 
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One important feature in the 7-e relation predicted by 
this model is that 7 ~ 1 for all < e < 1/6. For a pertur- 
bation with e in this range, the circular velocity of the halo 
increases rapidly with time, with a timescale that is shorter 
than a Hubble time. In this case, not only can particles with 
small apocenters (low orbital energies) reach the inner part 
of the halo, but also can many particles with large apoc- 
enters (high orbital energies) and small angular momenta. 
Velocity isotropisation mixes these orbits, resulting in 7 ~ 1 
as we have demonstrated. For e > 1/6 the gravitational po- 
tential is changing gradually and particles joining the halo 
have a similar energy and orbital shape J. The resulting pro- 
file can then be described by the self-similar solution that 
assumes the same orbital shape for all particles and 7 be- 
comes larger than one. Note that an inner logarithmic slope 
of approximately —1 results from a fast collapse and orbit 
isotropisation and that both conditions are required to pro- 
duce such an inner slope. If the collapse is fast (e ^ 1/6) but 
the velocity dispersion is not isotropic, the inner slope can 
be as shallow as (for constant J) and as steep as —2 (for 
radial infall). If the velocity dispersion is isotropic, but the 
mass accretion rate is small, i.e. e > 1, the inner slope can 
be much steeper than —1. 

The identification of this mechanism not only explains 
the approximate universality of the inner CDM halo slopes 
but also describes the variance observed in A''-body simula- 
tions. Gravitational collapse starts on small scales in a hi- 
erarchical model such as CDM. Deep dark-matter potential 
wells are created by subsequent non-linear gravitational col- 
lapse, but the potential well associated with a halo cannot 
deepen significantly during the slow accretion phase when 
the accretion time scale is longer than a Hubble time (Zhao 
et al. 2003a). Therefore, all halos must have gone through 
a phase of rapid accretion to establish their potential wells, 
even though different halos may have different mass accre- 
tion histories. Also, as shown in Li et al. (2005), the phase of 
rapid mass accretion is dominated by major mergers, which 
may effectively isotropise the velocity field. These are the 
two ingredients that are required to produce a p oc in- 
ner profile, and explain why such a inner profile results for 
halos with vastly different formation histories. If the poten- 
tial well of a halo is established at early times, the mass 
contained in the profile will be small, and the halo will 
have a high concentration since most of the mass accreted 
slowly. However, if a halo establishes its potential well re- 
cently, much of its mass will be in the profile, and the 
halo will have a low concentration. This correlation between 
halo concentration and the time of potential well formation 
matches cosmological TV-body simulations (e.g. Zhao et al. 
2003a;b). In extreme cases where the mass involved in the 
fast accretion is too small to be seen, we expect an inner 
profile steeper than r~^. 

So far we have only considered the origin of the inner 
profile. What about the outer r~'^ form in the universal 
profile? Consider a shell of mass M and of an initial radius 
ri. Suppose this mass shell collapses at a time t to a radius 
r. Assuming an Einstein-de Sitter Universe and neglecting 
the effect of shell crossing, we have r oc ri/Si{M) and t oc 
ti/6i{M)^^^ . Eliminating 5i{M) in these two relations and 
using M oc rf we obtain r oc M^''^f^/^. Thus, the density 
profile can be written as 



p(r) 



1 dM /dr 



47rr-2 dt 



M 



r» 2 + ^ 



(27) 



where = dlnM/dlnt. Now let us write the mass in two 
parts, M — Mg -\- AM, where Me is the mass of the halo at 
t — At, while AM is the mass accreted between time t — At 
and t. If AM increases as a power law of t, and if AM ^ Me, 
then p oc r~^. Since a mass shell settles into an equilibrium 
state over several dynamical times, the relevant scale for At 
is the dynamical time of the system. Thus, if there is a period 
of slow accretion, but where the total accreted mass is much 
smaller than the mass that has collapsed into the halo, an 



outer density profile with ■ 



results. However, the above 



argument also suggests that the outer density profile can 
change from halo to halo, depending on the mass accretion 
rate at the final stage of halo formation. For example, if M 
continues to grow as described by equation Q eventually a 
density profile with p oc in the outer parts will result. 

To demonstrate such a dependence of the outer profile 
on mass accretion history, we have carried out 1-D simu- 
lations for three cases with different recent mass accretion 
histories. These mass accretion histories are shown in the 
left panel of Figure |H1 and the resulting logarithmic slopes 
versus r are shown in the right panel. This figure shows 
that while halos with typical mass accretion histories have 
p oc outer density profiles, halos with slower accretion 
rates in the recent past (shown by the dashed and dotted 
curves) have steeper outer profiles. Therefore, the outer 
profile is not universal but a consequence of the form of mass 
accretion history of typical CDM halos and in the currently 
accepted cosmological model all halos will have a r~* outer 
profile in the distant future. 



5 SUMMARY AND DISCUSSION 

We have constructed a simple model that reproduces many 
of the properties of the CDM halo population. There are 
two essential ingredients: 1) a two-phase accretion history 
beginning with rapid accretion and potential-well deepening 
followed by slow accretion with little change in the poten- 
tial well; and 2) isotropisation of orbits during the phase 
of rapid growth. The density profiles obtained from various 
mass accretion histories all fit the NEW form. Our model 
also reproduces the correlation between the concentration 
and formation time observed in cosmological A^-body sim- 
ulations. In particular, the model predicts a roughly con- 
stant concentration, c ~ 5, for all halos that are still in 
the fast accretion phase, matching the results obtained by 
Zhao et al. (2003a;b). Combined with an ensemble of real- 
istic mass accretion histories parametrised from cosmologi- 
cal CDM simulations, the model reproduces the dependence 
of halo concentration on halo mass and the distribution of 
halo concentrations measured in these same cosmological N- 
body simulations. Our results demonstrate that the struc- 
tural properties of CDM halos are largely determined by 
their mass accretion histories. 

We can recover many of these results using a simple 
analytic model. Our model begins with scale-free perturba- 
tions of the form SM{r)/M{r) oc M(r)^'^ in an Einstein-de 
Sitter universe. Assuming an isotropic velocity dispersion, 
we find that the inner profile produced by the collapse of 
such perturbations is always a power law, p(r) oc r~^, with 
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Figure 8. The left panel shows three mass accretion histories that have different accretion rates in the recent past. The right panel 
shows, as a function of radius, the logarithmic slope of the density profile generated with each of these three mass accretion histories. 
Note that the outer slope can be as steep as —4 for a halo whose recent mass accretion rate is very small (dotted and dashed curves). 



7 = 1 for < e < 1/6, and 7 = 9e/(l + 3e) for e > 1/6. A 
model with < e < 1/6 has a rapidly deepening potential 
leading to orbit isotropy. This produces a shallower profile 
than a purely radial model. Assuming non-radial orbits of 
constant shape yields self-similar models but with a large 
range of possible inner slopes. The mixture resulting from 
isotropisation converges to a single inner slope p oc . This 
suggests that the inner profile of CDM halos is a nat- 
ural result of hierarchical models, where the potential well 
associated with a halo has to be built through a phase of 
rapid and violent accretion in which potential fiuctuations 
are expected to effectively isotropise the velocities of CDM 
particles. 

As mentioned in the introduction, a number of authors 
have previously noticed that tangential motions of parti- 
cles can cause flattening in the inner density profile of dark 
matter halos (e.g. White & Zartisky 1992; Ryden 1993; 
Sikivie et al. 1997; Subramanian et al. 2000; Subramanian 
2000; Hiotelis 2002; Le Delliou & Henriksen 2003; Shapiro 
et al. 2004; Barnes et al. 2005). Applying the collisionless 
Boltzmann equation to self-similar gravitational collapse, 
Subramanian (2000) and Subramanian et al. (2000) demon- 
strated that tangential velocities are required to obtain an 
inner profile that is shallower than that of a singular isother- 
mal sphere. They suggested that an appropriate mixture of 
radial and tangential velocities may produce an inner 
profile. However, these authors did not propose a specific 
model that predicts such a proflle. Furthermore, as we have 
shown in this paper, isotropic velocity dispersion alone is not 
suSicient for generating an inner profile; rapid accretion 
with a quickly deepening potential well is also required in 
our model. Thus, although our explanantion about the inner 
density profile of dark matter halos is related to those pre- 



sented in these earlier analyses, it provides additional phys- 
ical insight into the problem, in particular in connection to 
the properties of the mass accretion history that infiuence 
the final density profile. 

Since our model only relies on the mass accretion his- 
tory and orbit isotropisation to explain the origin of NFW 
profiles, it naturally explains why collapses with artificially 
reduced substructure (Moore et al. 1999) also lead to the 
universal profile, even when the initial condition is as smooth 
as three intersecting plane waves (Shapiro et al. 2004) . How- 
ever, if it is so smooth that the orbits are not isotropised then 
a steeper central slope may result. Note that although our 
explanation depends on different mass accretion histories it 
does not explicitly depend on the dynamical friction and 
tidal distribution of the substructures that makes up the 
actual mass accretion in cosmological A''-body simulations 
(e.g. Dekel et al. 2003). 

It is remarkable that our one-dimensional simulations 
with their two simple ingredients reproduce the structural 
features of the full three-dimensional simulations. In retro- 
spect, this result is consistent with the physical nature of 
hierarchical formation. Because dark-matter halos are ex- 
tended, even equal mass mergers are relatively quiescent, 
their mutual orbits slowly decaying by dynamical friction 
and ending with a low velocity merger, tidal dissolution, and 
phase mixing. Such events are likely to produce sufficient 
scattering to yield isotropisation but not as violent as envi- 
sioned by Lynden-Bell (1967). Therefore, once we have the 
two necessary ingredients, the final equilibrium profile may 
not be sensitive to the exact way in which the system settles 
into its final equilibrium configuration. Zhao et al. (2003a) 
show that a significant correlation still exists between the 
final and initial binding energies even for particles that are 
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accreted in the fast accretion phase based on 5 halos in their 
high-resolution simulations. Reinforcing this point, we also 
find that the final energy of a particle is correlated with its 
initial energy, and hence the energy distribution of particles 
is determined by the initial conditions rather than by com- 
plete relaxation. The match between our one-dimensional 
model and the three-dimensional cosmological simulations 
suggests that violent relaxation might not play an impor- 
tant role in redistributing the energies of particles except 
through isotropising the velocity field. 
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